Contents

Layout Layout

Using CAMS European air quality analysis from Copernicus Atmosphere Monitoring with RELIANCE services

Analysis over France

How to discover RELIANCE datacube resources (spatial and temporal search and subsetting)

This notebook shows how to discover and access the Copernicus Atmosphere Monitoring products available in the RELIANCE datacube resources, by using the functionalities provided in the Adam API . The process is structured in 6 steps, including example of data analysis and visualization with the Python libraries installed in the Jupyter environment.

import warnings

warnings.filterwarnings('ignore')

Step 1: Authentication

The following lines of code will show the personal Adam API-Key of the user and the endpoint currently in use, that provide access to the products in the related catalogue. At the end of the execution, if the authentication process is successfull the personal token and the expiration time should be returned as outputs.

!pip install adamapi

Requirement already satisfied: adamapi in /opt/conda/lib/python3.8/site-packages (2.0.11)
Requirement already satisfied: requests>=2.22.0 in /opt/conda/lib/python3.8/site-packages (from adamapi) (2.25.1)
Requirement already satisfied: imageio in /opt/conda/lib/python3.8/site-packages (from adamapi) (2.9.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (1.26.4)
Requirement already satisfied: idna<3,>=2.5 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (2.10)
Requirement already satisfied: certifi>=2017.4.17 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (2021.5.30)
Requirement already satisfied: chardet<5,>=3.0.2 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (4.0.0)
Requirement already satisfied: pillow in /opt/conda/lib/python3.8/site-packages (from imageio->adamapi) (8.1.2)
Requirement already satisfied: numpy in /opt/conda/lib/python3.8/site-packages (from imageio->adamapi) (1.20.2)

adam_key = open("adam-key").read().rstrip()

import adamapi as adam
a = adam.Auth()

a.setKey(adam_key)
a.setAdamCore('https://reliance.adamplatform.eu')
a.authorize() 

{'expires_at': '2021-10-14T10:18:21.777Z',
 'access_token': 'dae4a819e86f4b3086de49917de3c38b',
 'refresh_token': '98ff3616b6264bd683e7ed2ab5990ac8',
 'expires_in': 3600}

Step 2: Datasets Discovery

After the authorization, the user can browse the whole catalogue, structured as a JSON object after a pagination process, displaying all the available datasets. This operation can be executed with the getDatasets() function without including any argument. Some lines of code should be added to parse the Json object and extract the names of the datasets.The Json object can be handled as a Python dictionary.

Pre-filter datasets

We will discover all the available datasets in the ADAM platform but will only print elements of interest EU_CAMS e.g. European air quality datasets from Copernicus Atmosphere Monitoring Service

from adamapi import Datasets

datasets = Datasets(a)
catalogue = datasets.getDatasets()

#Extracting the size of the catalogue

total = catalogue['properties']['totalResults']
items = catalogue['properties']['itemsPerPage']
pages = total//items

print('----------------------------------------------------------------------')
print('\033[1m' + 'List of available datasets:')
print ('\033[0m')

#Extracting the list of datasets across the whole catalogue

for i in range(0,pages):
    page = datasets.getDatasets(page = i)
    for element in page['content']: 
        if 'EU_CAMS' in element['title'] :
            print(element['title'] + "\033[1m" + " --> datasetId "+ "\033[0m" + "= " + element['datasetId'])

----------------------------------------------------------------------
List of available datasets:

EU_CAMS_SURFACE_C2H3NO5_G --> datasetId = 69619:EU_CAMS_SURFACE_C2H3NO5_G
EU_CAMS_SURFACE_CO_G --> datasetId = 69620:EU_CAMS_SURFACE_CO_G
EU_CAMS_SURFACE_NH3_G --> datasetId = 69621:EU_CAMS_SURFACE_NH3_G
EU_CAMS_SURFACE_NMVOC_G --> datasetId = 69622:EU_CAMS_SURFACE_NMVOC_G
EU_CAMS_SURFACE_NO2_G --> datasetId = 69623:EU_CAMS_SURFACE_NO2_G
EU_CAMS_SURFACE_NO_G --> datasetId = 69624:EU_CAMS_SURFACE_NO_G
EU_CAMS_SURFACE_O3_G --> datasetId = 69625:EU_CAMS_SURFACE_O3_G
EU_CAMS_SURFACE_PM10_G --> datasetId = 69626:EU_CAMS_SURFACE_PM10_G
EU_CAMS_SURFACE_PM25_G --> datasetId = 69627:EU_CAMS_SURFACE_PM25_G
EU_CAMS_SURFACE_REC_G --> datasetId = 69628:EU_CAMS_SURFACE_REC_G
EU_CAMS_SURFACE_SIA_G --> datasetId = 69629:EU_CAMS_SURFACE_SIA_G
EU_CAMS_SURFACE_SO2_G --> datasetId = 69630:EU_CAMS_SURFACE_SO2_G
EU_CAMS_SURFACE_TEC_G --> datasetId = 69631:EU_CAMS_SURFACE_TEC_G

We are interested by Particulate matter < 10 µm (PM10) so we will discover EU_CAMS_SURFACE_PM10_G and print the metadata of this particular dataset, showing the data provenance.

datasetID = '69626:EU_CAMS_SURFACE_PM10_G'

print('\033[1;34m' + 'Metadata of ' + datasetID + ':')
print ('\033[0;0m')

paged = datasets.getDatasets(datasetID)
for i in paged.items():
    print("\033[1m" +  str(i[0]) + "\033[0m" + ': ' + str(i[1]))

Metadata of 69626:EU_CAMS_SURFACE_PM10_G:

datasetId: 69626:EU_CAMS_SURFACE_PM10_G
creationDate: 2021-07-12T02:00:00Z
dataType: Float32
epsg: 4326
keywords: []
license: {'documentationUrl': '', 'dataProviderName': 'ADS', 'dataProviderUrl': '', 'licenseId': '', 'dataPolicy': '', 'doi': '', 'credits': ''}
maxValue: [2645.995849609375]
minValue: [0.003789698239415884]
numberOfRecords: 57014
profile: {'profileSchema': 'eo_profile_schema.json', 'name': 'Earth Observation', 'mission': 'CAMS', 'sensor': 'CAMS', 'processingLevel': 'forecast', 'instrument': '', 'platform': ''}
resolutionUnit: degree
temporalResolution: Hourly
unit: 
unitDescription: 
updateDate: 2021-10-13T18:02:06Z
geometry: {'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}
resolutions: [0.1]
anyText: 4326,ESA,Daily,ESA
applications: ['Atmosphere']
datasetManager: govoni@meeo.it
datasetManagerOrganisation: meeo
description: This dataset provides daily air quality analyses and forecasts for Europe.

CAMS produces specific daily air quality analyses and forecasts for the European domain at significantly higher spatial resolution (0.1 degrees, approx. 10km) than is available from the global analyses and forecasts. The production is based on an ensemble of nine air quality forecasting systems across Europe. A median ensemble is calculated from individual outputs, since ensemble products yield on average better performance than the individual model products. The spread between the nine models are used to provide an estimate of the forecast uncertainty. The analysis combines model data with observations provided by the European Environment Agency (EEA) into a complete and consistent dataset using various data assimilation techniques depending upon the air-quality forecasting system used. In parallel, air quality forecasts are produced once a day for the next four days. Both the analysis and the forecast are available at hourly time steps at seven height levels.

Note that only nitrogen monoxide, nitrogen dioxide, sulphur dioxide, ozone, PM2.5, PM10 and dust are regularly validated against in situ observations, and therefore forecasts of all other variables are unvalidated and should be considered experimental.
title: EU_CAMS_SURFACE_PM10_G
noDataValue: 0
timeReferenceSystem: UTC
units: µg m-3
unitsDescription: PM10
services: ['MWCS', 'ADAM', 'AdamApi']
technicalManager: govoni@meeo.it
filtersEnabled: {'type': 'object', 'title': 'Dataset filters', 'properties': {'startDate': {'title': 'Start Date', 'type': 'string', 'format': 'date', 'text_rule': 'false', 'ops_only': 'false', 'math_rule': 'false', 'order_rule': 'false', 'pattern': '^d{4}-[01]d-[0-3]d(T[0-2]d:[0-9]d:[0-9]dZ)?$'}, 'endDate': {'title': 'End Date', 'type': 'string', 'format': 'date', 'ops_only': 'false', 'text_rule': 'false', 'math_rule': 'false', 'order_rule': 'false', 'pattern': '^d{4}-[01]d-[0-3]d(T[0-2]d:[0-9]d:[0-9]dZ)?$'}}}
starDate: 2018-07-12T00:00:00Z
endDate: 2021-10-12T23:00:00Z
accounting: {'unit': 0.0032}

Step 3: Products Discovery

The products discovery operation related to a specific dataset is implemented in the Adam API with the getProducts() function. A combined spatial and temporal search can be requested by specifying the datasetId for the selected dataset,the geometry argument that specifies the Area Of Interest and a temporal range, defined by startDate and endDate . The geometry must always be defined by a GeoJson object that describes the polygon in the counterclockwise winding order. The optional arguments startIndex and maxRecords can set the list of the results returned as an output. The results of the search are displayed with their metadata and they are sorted starting from the most recent product.

Generate geometry for France

The geometry field is extracted from a GeoJSON object , retrieving the value of the “feature” element.

!wget https://raw.githubusercontent.com/mledoze/countries/master/data/fra.geo.json

--2021-10-14 09:18:26--  https://raw.githubusercontent.com/mledoze/countries/master/data/fra.geo.json
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.109.133, 185.199.108.133, 185.199.111.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.109.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 42936 (42K) [text/plain]
Saving to: ‘fra.geo.json.2’

fra.geo.json.2      100%[===================>]  41.93K  --.-KB/s    in 0.001s  

2021-10-14 09:18:26 (46.5 MB/s) - ‘fra.geo.json.2’ saved [42936/42936]

Search data

!pip install geojson_rewind

Requirement already satisfied: geojson_rewind in /opt/conda/lib/python3.8/site-packages (1.0.1)

from adamapi import Search
from geojson_rewind import rewind
import json

The GeoJson object needs to be rearranged according to the counterclockwise winding order.This operation is executed in the next few lines to obtain a geometry that meets the requirements of the method. Geom_1 is the final result to be used in the discovery operation.

with open('fra.geo.json') as f:
    geom_dict = json.load(f)
output = rewind(geom_dict)    
geom_1 = str(geom_dict['features'][0]['geometry'])

Copernicus air quality analyses are hourly product but when we select a given date, we will only get the first 10 products. Below, we make a list of the first 10 available products for the 1st September 2020 e.g. we restrict our search to this date.

start_date = '2019-09-01'
end_date = '2019-09-01'

search = Search( a )
results = search.getProducts(
    datasetID, 
    geometry= geom_1,
    startDate = start_date,
    endDate = end_date
 )

# Printing the results

print('\033[1m' + 'List of available products:')
print ('\033[0m')
count = 1
for i in results['content']:

        print("\033[1;31;1m" + "#" + str(count))
        print ('\033[0m')
        for k in i.items():
            print(str(k[0]) + ': ' + str(k[1]))
        count = count+1
        print('------------------------------------')

List of available products:

#1

_id: {'$oid': '60ef2c3d6dfebc0806c75279'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T23:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901230000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#2

_id: {'$oid': '60ef2c3d6dfebc0806c75276'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T22:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901220000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#3

_id: {'$oid': '60ef2c3d6dfebc0806c7526b'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T21:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901210000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#4

_id: {'$oid': '60ef2c3d6dfebc0806c75269'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T20:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901200000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#5

_id: {'$oid': '60ef2c3d6dfebc0806c75260'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T19:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901190000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#6

_id: {'$oid': '60ef2c3d6dfebc0806c75263'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T18:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901180000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#7

_id: {'$oid': '60ef2c3d6dfebc0806c75255'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T17:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901170000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#8

_id: {'$oid': '60ef2c3d6dfebc0806c75250'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T16:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901160000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#9

_id: {'$oid': '60ef2c3d6dfebc0806c7524a'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T15:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901150000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------
#10

_id: {'$oid': '60ef2c3d6dfebc0806c7524c'}
datasetId: EU_CAMS_SURFACE_PM10_G
productDate: 2019-09-01T14:00:00Z
productId: EU_CAMS_SURFACE_PM10_G_20190901140000.tif
geometry: {'geometries': [{'type': 'Polygon', 'coordinates': [[[-25.000012, 29.999997], [44.999988, 29.999997], [44.999988, 71.999997], [-25.000012, 71.999997], [-25.000012, 29.999997]]]}], 'type': 'GeometryCollection'}
insertDate: 2021-07-14T18:26:05Z
status: Online
------------------------------------

Step 4: Data Access

After the data discovery operation that retrieves the availability of products in the catalogue, it is possible to access the data with the getData function. Each product in the output list intersects the selected geometry and the following example shows how to access a specific product from the list of results obtained in the previous step. While the datasetId is always a mandatory parameter, for each data access request the getData function needs only one of the following arguments: geometry or productId , that is the value of the _id field in each product metadata. In the case of a spatial and temporal search the geometry must be provided to the function, together with the time range of interest. The output of the getData function is always a .zip file containing the data retrieved with the data access request, providing the spatial subset of the product. The zip file will contain a geotiff file for each of the spatial subsets extracted in the selected time range.

from adamapi import GetData

Define a function to select a time range and get data

import pathlib

def getZipData(auth, dataset_info):
    if not (pathlib.Path(pathlib.Path(dataset_info['outputFname']).stem).exists() or pathlib.Path(dataset_info['outputFname']).exists()):
        data=GetData(auth)
        image = data.getData(
        datasetId = dataset_info['datasetID'],
        startDate = dataset_info['startDate'],
        endDate = dataset_info['endDate'],
        geometry = dataset_info['geometry'],
        outputFname = dataset_info['outputFname'])
        print(image)

Get PM10 for each day of September 2019, 2020 and 2021 (time 00:00:00)

This process can take a bit of time so be patient!

import time
from IPython.display import clear_output

start = time.time()

for year in ['2019', '2020', '2021']:
    datasetInfo = {
    'datasetID' : '69626:EU_CAMS_SURFACE_PM10_G',
    'startDate' : year + '-09-01',
    'endDate' : year + '-09-30',
    'geometry' : geom_1,
    'outputFname' : './PM10_FR_ADAMAPI_' + year + '.zip'
    }
    getZipData(a, datasetInfo)
    
end = time.time()
clear_output(wait=True)
delta1 = end - start
print('\033[1m'+'Processing time: ' + str(round(delta1,2)))

Processing time: 0.0

Step 5: Data Analysis and Visualization

The data retrieved via the Adam API are now available as a zip file that must be unzipped to directly handle the data in a geotiff format. Then with the Python packages provided in the Jupyter environment it is possible to process and visualized the requested product.

Unzip data

import zipfile

def unzipData(filename):
    with zipfile.ZipFile(filename, 'r') as zip_ref:
        zip_ref.extractall(path = pathlib.Path(filename).stem)

for year in ['2019', '2020', '2021']:
    filename = './PM10_FR_ADAMAPI_' + year + '.zip'
    if not pathlib.Path(pathlib.Path(filename).stem).exists():
        unzipData(filename)

Read data and make a monthly average

import xarray as xr
import glob

We now read these files using xarray. First, we make a list of all the geotiff files in a given folder. To ensure each raster is labelled correctly with its time, we can use a helper function paths_to_datetimeindex() to extract time information from the file paths we obtained above. We then load and concatenate each dataset along the time dimension using xarray.open_rasterio(), convert the resulting xarray.DataArray to a xarray.Dataset, and give the variable a more useful name (PM10)

from datetime import datetime

def paths_to_datetimeindex(paths):
    return  [datetime.strptime(date.split('_')[-1].split('.')[0], '%Y-%m-%dt%f') for date in paths]

def getData(dirtif):
    geotiff_list = glob.glob(dirtif)
    # Create variable used for time axis
    time_var = xr.Variable('time', paths_to_datetimeindex(geotiff_list))
    # Load in and concatenate all individual GeoTIFFs
    geotiffs_da = xr.concat([xr.open_rasterio(i, parse_coordinates=True) for i in geotiff_list],
                        dim=time_var)
    # Covert our xarray.DataArray into a xarray.Dataset
    geotiffs_da = geotiffs_da.to_dataset('band')
    # Rename the dimensions to make it CF-convention compliant
    geotiffs_da = geotiffs_da.rename_dims({'y': 'latitude', 'x':'longitude'})
    # Rename the variable to a more useful name
    geotiffs_da = geotiffs_da.rename_vars({1: 'PM10', 'y':'latitude', 'x':'longitude'})

    return geotiffs_da

geotiff_ds = getData('./PM10_FR_ADAMAPI_20*/*.tif')
geotiff_ds['PM10'].attrs = {'units' : 'µg m-3', 'long_name' : 'Particulate matter < 10 µm'}
geotiff_ds

<xarray.Dataset>
Dimensions:    (time: 90, latitude: 98, longitude: 148)
Coordinates:
  * latitude   (latitude) float64 51.05 50.95 50.85 50.75 ... 41.55 41.45 41.35
  * longitude  (longitude) float64 -5.15 -5.05 -4.95 -4.85 ... 9.35 9.45 9.55
  * time       (time) datetime64[ns] 2020-09-19 2020-09-20 ... 2021-09-12
Data variables:
    PM10       (time, latitude, longitude) float32 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    transform:      (0.1, 0.0, -5.200012009298998, 0.0, -0.1, 51.09999699761336)
    crs:            +init=epsg:4326
    res:            (0.1, 0.1)
    is_tiled:       0
    nodatavals:     (0.0,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area

Analyze data

Make yearly average for September month

geotiff_dm = geotiff_ds.groupby('time.year').mean('time', keep_attrs=True)

geotiff_dm

<xarray.Dataset>
Dimensions:    (year: 3, latitude: 98, longitude: 148)
Coordinates:
  * latitude   (latitude) float64 51.05 50.95 50.85 50.75 ... 41.55 41.45 41.35
  * longitude  (longitude) float64 -5.15 -5.05 -4.95 -4.85 ... 9.35 9.45 9.55
  * year       (year) int64 2019 2020 2021
Data variables:
    PM10       (year, latitude, longitude) float32 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    transform:      (0.1, 0.0, -5.200012009298998, 0.0, -0.1, 51.09999699761336)
    crs:            +init=epsg:4326
    res:            (0.1, 0.1)
    is_tiled:       0
    nodatavals:     (0.0,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area

Visualize data

!pip install cmaps

Requirement already satisfied: cmaps in /opt/conda/lib/python3.8/site-packages (1.0.4)
Requirement already satisfied: matplotlib in /opt/conda/lib/python3.8/site-packages (from cmaps) (3.4.3)
Requirement already satisfied: numpy in /opt/conda/lib/python3.8/site-packages (from cmaps) (1.20.2)
Requirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (2.8.1)
Requirement already satisfied: pyparsing>=2.2.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (2.4.7)
Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (0.10.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (1.3.1)
Requirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (8.1.2)
Requirement already satisfied: six in /opt/conda/lib/python3.8/site-packages (from cycler>=0.10->matplotlib->cmaps) (1.16.0)

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cmaps

# generate figure
proj_plot = ccrs.Mercator(central_longitude=24.0)

lcmap = cmaps.BlueYellowRed
# Only plot values greater than 0
p = geotiff_dm['PM10'].where(geotiff_dm['PM10'] > 0).plot(x='longitude', y='latitude',transform=ccrs.PlateCarree(),
                            subplot_kws={"projection": proj_plot},
                            size=8,
                            col='year', col_wrap=3, robust=True, cmap=lcmap, add_colorbar=True)

# We have to set the map's options on all four axes
for ax,i in zip(p.axes.flat,  geotiff_dm.year.values):
    ax.coastlines()
    ax.set_title('Surface PM10\n' + 'September ' + str(i), fontsize=10)

plt.savefig('b2drop/PM10_september_FR_2019-2021.png')
../_images/RELIANCE_v0.2-France_58_0.png

Plot one single date

fig=plt.figure(figsize=(10,10))
# Define the projection
crs=ccrs.PlateCarree()

# We're using cartopy and are plotting in Orthographic projection 
# (see documentation on cartopy)
ax = plt.subplot(1, 1, 1, projection=ccrs.Mercator(central_longitude=12.0))
ax.coastlines(resolution='10m')

# custom colormap

lcmap = cmaps.BlueYellowRed

# We need to project our data to the new Mercator projection and for this we use `transform`.
# we set the original data projection in transform (here PlateCarree)
# we only plot values greather than 0
img = geotiff_ds['PM10'].where(geotiff_ds['PM10'] > 0).sel(time='2021-09-15').plot(ax=ax, transform=ccrs.PlateCarree(), cmap=lcmap)  

# Title for plot
plt.title('Surface PM10 \n 15th September 2021 over France',fontsize = 16, fontweight = 'bold', pad=10)
plt.savefig('b2drop/PM10_september_FR_2021-09-15.png')
../_images/RELIANCE_v0.2-France_60_0.png

Interactive plot with bokeh

import hvplot.xarray

PM10 over France from 1st September to 30th September 2021

geotiff_ds = geotiff_ds.sortby('time')

geotiff_ds.hvplot.labels( text_font_size='6pt', text_color='blue').opts(xoffset=20)
geotiff_ds.where(geotiff_ds['PM10'] > 0).sel(time=slice('2021-09-01', '2021-09-30')).PM10.hvplot(
    groupby="time",
    cmap=lcmap,
    geo=True,
    coastline='10m',
    frame_width=400,
    clim=(0,35) 
)

Save Data Cube selection into netCDF

geotiff_ds.to_netcdf("b2drop/PM10_september_FR_2019-2021.nc")

Plot a timeseries

PM10 over Paris from 1st September to 30th September 2021

Sort times and then select one location

latitude = 48.864716
longitude = 2.349014

ts = geotiff_ds['PM10'].where(geotiff_ds['PM10'] > 0).sel(latitude=latitude, longitude=longitude, method='nearest').sel(time=slice('2021-09-01', '2021-09-30'))

ts.hvplot(color='purple') * ts.hvplot.scatter(marker='o', color='blue', size=15)

Save timeseries into tabular (csv) file

ts.to_series().to_csv("b2drop/timeseries_paris_september2021.csv", sep='\t')

Create Research Object in ROHUB

import os
import pathlib
from modules.rohub import rohub, settings

Authenticating

  • If the code cell below fails, make sure you have created the two files:

    • rohub-user: contains your rohub username

    • rohub-pwd: contains your rohub password

rohub_user = open("rohub-user").read().rstrip()
rohub_pwd = open("rohub-pwd").read().rstrip()

rohub.login(username=rohub_user, password=rohub_pwd)

Logged successfully as annefou.

Create new basic RO

ro_title="Copernicus Atmosphere Monitoring Service Code Research Object"
ro_research_areas=["Earth sciences"]

ro = rohub.ros_create(title=ro_title, research_areas=ro_research_areas)

Research Object was successfully created with id = 930fcc3c-be00-4a17-a395-d9c1cba43cd5

Show basic metadata

ro.show_metadata()

{'identifier': '930fcc3c-be00-4a17-a395-d9c1cba43cd5',
 'shared_link': 'https://w3id.org/ro-id/930fcc3c-be00-4a17-a395-d9c1cba43cd5',
 'title': 'Copernicus Atmosphere Monitoring Service Code Research Object',
 'description': None,
 'status': 'LIVE',
 'access_mode': 'PUBLIC',
 'type': 'Basic Research Object',
 'template': None,
 'created': '2021-10-14T09:18:58.335017Z',
 'creator': 'annefou',
 'modificator': 'annefou',
 'modified': None,
 'importer': None,
 'research_areas': ['Earth sciences'],
 'owner': 'annefou',
 'editors': [],
 'readers': [],
 'rating': '0.00',
 'number_of_ratings': 0,
 'number_of_likes': 0,
 'number_of_dislikes': 0,
 'quality': 0,
 'size': 0,
 'creation_mode': 'MANUAL',
 'doi': None,
 'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/ros/930fcc3c-be00-4a17-a395-d9c1cba43cd5/'}

Add metadata associated to our Research Object

ro.description="This Research Object demonstrate how to use CAMS European air quality analysis from Copernicus Atmosphere Monitoring with RELIANCE services and compute monthly map of PM10 over a given geographical area."
ro.ros_type="Code Research Research Object"

ro.update()

{'identifier': '930fcc3c-be00-4a17-a395-d9c1cba43cd5',
 'shared_link': 'https://w3id.org/ro-id/930fcc3c-be00-4a17-a395-d9c1cba43cd5',
 'title': 'Copernicus Atmosphere Monitoring Service Code Research Object',
 'description': 'This Research Object demonstrate how to use CAMS European air quality analysis from Copernicus Atmosphere Monitoring with RELIANCE services and compute monthly map of PM10 over a given geographical area.',
 'status': 'LIVE',
 'access_mode': 'PUBLIC',
 'type': 'Code Research Research Object',
 'template': None,
 'created': '2021-10-14T09:18:58.335017Z',
 'creator': 'annefou',
 'modificator': 'annefou',
 'modified': '2021-10-14T09:19:01.960692Z',
 'importer': None,
 'research_areas': ['Earth sciences'],
 'owner': 'annefou',
 'editors': [],
 'readers': [],
 'rating': '0.00',
 'number_of_ratings': 0,
 'number_of_likes': 0,
 'number_of_dislikes': 0,
 'quality': 0,
 'size': 1732,
 'creation_mode': 'MANUAL',
 'doi': None,
 'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/ros/930fcc3c-be00-4a17-a395-d9c1cba43cd5/'}

Aggregate resource (internal)

Add our Jupyter Notebook to our Research Object

resi_file_path=os.path.join(os.environ['HOME'], 'RELIANCE_v0.2-France.ipynb')

if pathlib.Path(resi_file_path).exists():
    resi_res_type="Python Script"
    resi_title="Compute and visualize monthly mean maps of CAMS PM10 for different years"
    resi_description="Jupyter Notebook for discovering, accessing and processing and visualizing RELIANCE data cube of CAMS Particule matter < 10 μm [μg/m3]"
    ro.add_internal_resource(res_type=resi_res_type,file_path=resi_file_path, title=resi_title, description=resi_description)

Aggregate resource (external)

Add Results to our Research Object

  • We will be adding all the resources generated by our notebook (data and plots)

  • All our data and plots are shared in B2DROP so we will get the shared link from B2DROP and add it to our research object

Add first plot as sketch to my Research Object (external resource from B2DROP)

rese_file_url = getB2dropLink('PM10_september_FR_2019-2021.png', b2drop_user, b2drop_pwd)
rese_res_type="Sketch"
rese_title="Particule matter < 10 μm [μg/m3] over France for September 2019, 2020 and 2021"
rese_description="Monthly average maps of CAMS Particule matter < 10 μm [μg/m3] over France in 2019, 2020 and 2021"
ro.add_external_resource(res_type=rese_res_type, url=rese_file_url, title=rese_title, description=rese_description)

{'identifier': '5d46557d-43da-46ca-bcc4-18f9a3fa364f',
 'title': 'Particule matter < 10 μm [μg/m3] over France for September 2019, 2020 and 2021',
 'folder': None,
 'ros': ['930fcc3c-be00-4a17-a395-d9c1cba43cd5'],
 'description': 'Monthly average maps of CAMS Particule matter < 10 μm [μg/m3] over France in 2019, 2020 and 2021',
 'url': 'https://b2drop.eudat.eu/s/HJDimSrpjqTDCxQ/preview',
 'name': 'Particule matter < 10 μm [μg/m3] over France for September 2019, 2020 and 2021',
 'filename': None,
 'path': 'https://b2drop.eudat.eu/s/HJDimSrpjqTDCxQ/preview',
 'size': None,
 'download_url': 'https://b2drop.eudat.eu/s/HJDimSrpjqTDCxQ/preview',
 'type': 'Sketch',
 'doi': None,
 'read_only': False,
 'cloned': None,
 'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/resources/5d46557d-43da-46ca-bcc4-18f9a3fa364f/'}

Add second plot as Image to my Research Object (external resource from B2DROP)

rese_file_url = getB2dropLink('PM10_september_FR_2021-09-15.png', b2drop_user, b2drop_pwd)
rese_res_type="Image"
rese_title="Particule matter < 10 μm [μg/m3] over France on September 15, 2021"
rese_description="Daily average maps of CAMS Particule matter < 10 μm [μg/m3] over France on September 15, 2021"
ro.add_external_resource(res_type=rese_res_type, url=rese_file_url, title=rese_title, description=rese_description)

{'identifier': '4ea9eeea-1924-4e08-8c21-d249544a5d34',
 'title': 'Particule matter < 10 μm [μg/m3] over France on September 15, 2021',
 'folder': None,
 'ros': ['930fcc3c-be00-4a17-a395-d9c1cba43cd5'],
 'description': 'Daily average maps of CAMS Particule matter < 10 μm [μg/m3] over France on September 15, 2021',
 'url': 'https://b2drop.eudat.eu/s/StZTWgZpBYd3qtR/preview',
 'name': 'Particule matter < 10 μm [μg/m3] over France on September 15, 2021',
 'filename': None,
 'path': 'https://b2drop.eudat.eu/s/StZTWgZpBYd3qtR/preview',
 'size': None,
 'download_url': 'https://b2drop.eudat.eu/s/StZTWgZpBYd3qtR/preview',
 'type': 'Image',
 'doi': None,
 'read_only': False,
 'cloned': None,
 'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/resources/4ea9eeea-1924-4e08-8c21-d249544a5d34/'}

Add timeseries for Paris (external resource from B2DROP)

rese_file_url = getB2dropLink('timeseries_paris_september2021.csv', b2drop_user, b2drop_pwd)
rese_res_type="Result"
rese_title="Timeseries of particule matter < 10 μm [μg/m3] over Paris in september 2021"
rese_description="Daily average of CAMS Particule matter < 10 μm [μg/m3] over Paris in September 2021"
ro.add_external_resource(res_type=rese_res_type, url=rese_file_url, title=rese_title, description=rese_description)

{'identifier': '0f1bd850-2cfb-4754-ae7d-b7f1e58a2b66',
 'title': 'Timeseries of particule matter < 10 μm [μg/m3] over Paris in september 2021',
 'folder': None,
 'ros': ['930fcc3c-be00-4a17-a395-d9c1cba43cd5'],
 'description': 'Daily average of CAMS Particule matter < 10 μm [μg/m3] over Paris in September 2021',
 'url': 'https://b2drop.eudat.eu/s/YrC9GWZyrjK4GHr/preview',
 'name': 'Timeseries of particule matter < 10 μm [μg/m3] over Paris in september 2021',
 'filename': None,
 'path': 'https://b2drop.eudat.eu/s/YrC9GWZyrjK4GHr/preview',
 'size': None,
 'download_url': 'https://b2drop.eudat.eu/s/YrC9GWZyrjK4GHr/preview',
 'type': 'Result',
 'doi': None,
 'read_only': False,
 'cloned': None,
 'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/resources/0f1bd850-2cfb-4754-ae7d-b7f1e58a2b66/'}

Add netCDF file corresponding to Data cube selection (external resource from B2DROP)

rese_file_url = getB2dropLink('PM10_september_FR_2019-2021.nc', b2drop_user, b2drop_pwd)
rese_res_type="Result"
rese_title="netCDF data for daily PM10 concentration over France in September 2019, 2020 and 2021"
rese_description="netCDF data corresponding to daily average of CAMS Particule matter < 10 μm [μg/m3] over France for September 2019, September 2020 and September 2021"
ro.add_external_resource(res_type=rese_res_type, url=rese_file_url, title=rese_title, description=rese_description)

{'identifier': '07925723-f747-4f2e-aca6-d9aec1a36885',
 'title': 'netCDF data for daily PM10 concentration over France in September 2019, 2020 and 2021',
 'folder': None,
 'ros': ['930fcc3c-be00-4a17-a395-d9c1cba43cd5'],
 'description': 'netCDF data corresponding to daily average of CAMS Particule matter < 10 μm [μg/m3] over France for September 2019, September 2020 and September 2021',
 'url': 'https://b2drop.eudat.eu/s/Bp4J7BgEgGA3KAi/preview',
 'name': 'netCDF data for daily PM10 concentration over France in September 2019, 2020 and 2021',
 'filename': None,
 'path': 'https://b2drop.eudat.eu/s/Bp4J7BgEgGA3KAi/preview',
 'size': None,
 'download_url': 'https://b2drop.eudat.eu/s/Bp4J7BgEgGA3KAi/preview',
 'type': 'Result',
 'doi': None,
 'read_only': False,
 'cloned': None,
 'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/resources/07925723-f747-4f2e-aca6-d9aec1a36885/'}

Add geolocation to my Research Object

geolocation_file_path = os.path.join(os.environ['HOME'],  'fra.geo-ld.json') 
ro.add_geolocation(body_specification_json=geolocation_file_path)

[{'identifier': 'b2efe8e0-5706-423c-949a-2bebb719998c',
  'name': 'annotations/b2efe8e0-5706-423c-949a-2bebb719998c',
  'filename': '.ro/annotations/b2efe8e0-5706-423c-949a-2bebb719998c.ttl',
  'relation': False,
  'read_only': False,
  'target': '930fcc3c-be00-4a17-a395-d9c1cba43cd5',
  'ro': '930fcc3c-be00-4a17-a395-d9c1cba43cd5',
  'resources': [],
  'created': '2021-10-14T09:19:14.586237Z',
  'creator': 'annefou',
  'modificator': None,
  'modified': None,
  'api_link': 'https://rohub2020-rohub.apps.paas-dev.psnc.pl/api/annotations/b2efe8e0-5706-423c-949a-2bebb719998c/'}]

Publish my Research Object to Zenodo

  • Take a snapshot and publis my RO to zenodo

snapshot_id=ro.snapshot(publication_services=["Zenodo"])

Task done successfully

Load the published Research Object

published_ro = rohub.ros_load(identifier=snapshot_id)

Research Object was successfully loaded with id = 2aec1040-82ee-4107-9394-d72a24c9834a

previous

Using CAMS European air quality analysis from Copernicus Atmosphere Monitoring with RELIANCE services

next

Notebooks from researchers